function Q = portmanteauTestGarch(x,nlag)

    n = length(x);    
%     h = abs(h);
    x = x - mean(x);
    

%     mycov = @(h,k)(sum((x(1:n-k).^2) .*x(k+1:end).*x(k+1-h:n-h))
    % k > h
%     mycov = @(h,k)(sum( (x(max(h,k)+1 : n).^2) .* x(1 : n-max(h,k) ) ...
%                          .*x(max(h,k)+1-min(h,k) : n-min(h,k))) );
%     covMatrix = arrayfun(@(h,k) mycov2(x,h,k),(1:nlag), (1:nlag),'UniformOutput',true);
    Q = zeros(nlag);    
    for h = 1:nlag
        for k = h : nlag
            Q(h,k) = mycov2(x,h,k);        
        end
    end    
end

function xx = mycov2(x,h,k)
% calculate the covariance(i,j): (x_{t})^2 .* x_{t-k} .* x_{t-h}
    if(h == k)
        fprintf('..');
    end
    n = length(x);
    x1 = (x(max(h,k)+1 : n).^2);
    x2 = x(1 : n-max(h,k) );
    x3 = x( max(h,k) + 1 - min(h,k) : n - min(h,k));
    
    xx = sum(x1.*x2.*x3)./n;
   
end